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Abstract 


A detailed description of the development of the tangent linear model (TLM) and its 
adjoint model of the Relaxed Arakawa-Schubert moisture parameterization package used 
in NASA GEOS-1 C-Grid GCM (Version 5.2) is presented. The notational conventions 
used in the TLM and its adjoint codes are described in detail. 



Contents 


List of Figures vii 

1 Introduction 1 

2 Description of the Moisture Parameterization Scheme 1 

2.1 Cumulus Parameterization Package of the NASA GEOS-1 GCM 5 

2.2 Description of the Discretization of the Moisture Physics Parameterization . 12 

3 Tangent Linear Model of the Moist Process Physics Package 17 

3.1 Linearized Discrete Dynamical Equations 23 

3.2 Coding of the Tangent Linear Model 23 

3.3 Notational Convention for Variables and Subroutines Used in the Tangent 

Linear Model Code 24 

4 Adjoint Model of the Moist Process Physics Package 24 

4.1 Using the Adjoint Method to Calculate the Gradient of a Cost Function . . 24 

4.2 Coding of the Adjoint Model 25 

4.3 Notational Convention for Variables and Subroutines Used in the Adjoint 

Model Code 26 

4.4 Verification of the Correctness of the Adjoint Model 26 

Acknowledgments 31 

References 33 


v 



List of Figures 


1 Variation of the (f>(a) with respect to logo (gradient check of correctness of 

adjoint model). Integration period is 6 hours and t = 6 hours model generated 
observations were used. January 1, 1985 00Z DAO’s data was used as t = 0 
observations. The first guess is the shifted 6-hour initial condition. The time 
integration scheme employed is the leapfrog scheme 28 

2 As in Figure 1, but for the Matsuno time integration scheme 28 

3 Variation of the log|<j!>(o) — 1| with respect to logo. Integration period is 6 
hours and t = 6 hours model generated observations were used. January 1, 

1985 00Z DAO’s data is used as t = 0 observations. The first guess is the 
shifted 6-hour initial condition. The time integration scheme employed is the 
leapfrog scheme 29 

4 As in Figure 3, but for the Matsuno time integration scheme 29 


vii 



1 Introduction 


The GEOS-1 C-Grid GCM was developed by the Data Assimilation Office (DAO) at God- 
dard Laboratory for Atmosphere (GLA), NASA/GSFC (Takacs, et ah, 1994; Suarez and 
Takacs 1994) to be used in conjunction with an analysis scheme to produce a multi-year 
global atmospheric data set for climate research (Schubert et ah, 1993). It has an advanced 
structure, i.e., a “plug-compatible” structure. It means that if “plug-compatible” rules are 
followed in coding different GCMs and parameterizations, codes can be “unplugged” from 
one model and “plugged” into another with little coding effort. Thus each part of the 
GEOS-1 C-grid GCM can be used independently at another GCM. For instance, the full 
physics package of GEOS-1 C-grid GCM has been used into NASA/GLA Semi-Lagrangian 
Semi-Implicit (SLSI) GCM. The DAO and the Climate and Radiation Branch at GLA, 
NASA/GSFC have produced a library of physical parameterizations and dynamical algo- 
rithms which may be utilized for various GCM applications. 

There are four physics parameterization packages in the GEOS-1 GCM, i.e., the Relaxed 
Arakawa-Schubert (RAS) subgrid moisture parameterization and large-scale convection 
schemes, the shortwave radiation scheme, the longwave radiation scheme, and the tur- 
bulence parameterization scheme. Recently, the land surface model (LSM) described by 
Koster and Suarez (1992) has been fully coupled to the Aries GCM at NASA/GSFC, and 
the coupled models have been used to address a number of climate-related problems (Koster 
and Suarez 1996). In the physics packages, the moisture process plays an essential role to- 
wards improving the quality of the GCM productions. 

In order to obtain a 4-D variational data assimilation system based on the NASA GEOS-1 C- 
grid GCM including the moisture parameterization process, we need to develop the tangent 
linear and corresponding adjoint model of this physics package. This document describes the 
development of the tangent linear and adjoint models of the RAS and large-scale convection 
schemes as well as the formation of clouds which provides the cloud information used for 
cloud-radiative interactions. 

In Section 2 we provide a condensed description of the moisture parameterization package, 
including the basic original equations and their discrete forms. Section 3 describes and 
documents in detail the derivation and coding at the tangent linear model (TLM), and the 
notational conventions used. Section 4 describes the derivation of the adjoint model code 
and discusses its validity. 


2 Description of the Moisture Parameterization Scheme 


In a large-scale disturbance, there are many individual cumulus clouds whose spatial and 
temporal scales are much smaller than the disturbance itself. Because of this scale separa- 
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tion, it may be possible to predict the time evolution with the collective influence of all the 
smaller scale cumulus clouds. This is the goal of cumulus parameterization. 

Subgrid scale cumulus convection parameterization is a crucial component of any GCMs. 
Without including it in a GCM, one cannot hope to simulate the right general circula- 
tion patterns especially over the tropical region, and some other regions, where cumulus 
convective activities are frequent and strong. 

The necessity of parameterization of cumulus convection was understood and became clear 
after early research efforts failed to explain theoretically the size and growth rate of tropical 
cyclones (Lilly 1960; Yanai 1964). Among of the early parameterization researches, Charney 
and Elisaaen (1964) and Ooyama (1964) presented their classical papers in which the concept 
of “conditional instability of the second kind” (CISK) was first introduced. CISK describes 
the cooperative interaction between the cumulus scale and large scale, i.e., the large-scale 
circulation is responsible for organizing and maintaining cumulus convection by providing 
the necessary horizontal transport of water vapor, while the cumulus-scale drives the large- 
scale circulation through the release of latent heat in deep convective elements. CISK 
mechanism treats cumulus activities to be a function of the large-scale fields and it is now 
commonly referred to as cumulus parameterization. The simple parameterizations used in 
these papers led to considerable success in the numerical simulation of tropical cyclones 
(e.g., Ooyama 1969). Due to the high degree of empiricism and intuition, and the lack of 
theoretical framework for describing the mutual interaction between a cumulus ensemble 
and the large-scale environment, these early parameterizations were too crude to be used 
in GCM. 

Ooyama (1971) first proposed a cumulus parameterization theory which took into account 
the coexistence of a spectrum of clouds. He assumed that cumulus clouds can be repre- 
sented as non-interacting spherical bubbles, dispatched from the mixed layer. He concluded 
that the problem of parameterization of cumulus convection reduces to a determination of 
the dispatch function thus his parameterization scheme is not closed due to the unknown 
dispatch function. 

Before 1974, the basic physical image related to the cumulus parameterization, was that 
an existing cumulus cloud ensemble produces time changes in the large-scale temperature 
and moisture fields (Araka.wa. 1969, 1971, 1972; Betts 1973a, b; Gray 1972; Lopez 1972a, b; 
Ooyama. 1971; Yanai 1971a., b; Yanai et al. 1973). Cumulus convection modifies the large- 
scale temperature and moisture fields through detrainment and cumulus-induced subsidence 
in the environment. The detrainment causes large-scale cooling and moistening, and the 
cumulus-induced subsidence causes large-scale warming and drying. Based on 1956 Marshall 
Island observation data., Yanai et al. (1973) quantitatively derived these effects from a 
combination of observed large-scale heat and moisture budget over an area covered by 
cloud cluster and a. cumulus ensemble model which is similar to the one used in Arakawa 
and Schubert (1974). Their results show the importance of coexistence of shallow clouds 
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with deep clouds in maintaining the large-scale heat and moisture budgets. 

Based on the research results of the cumulus parameterization mentioned above, Arakawa 
and Schubert (1974, hereafter AS) presented a remarkable achievement, namely they pro- 
posed a closed cumulus parameterization theory that describes the mutual interaction of 
an ensemble of cumulus clouds with the large-scale environment. This is perhaps the most 
physically complete approach to the issue of cumulus parameterization up to now. The 
cloud ensemble is represented by a spectrum of idealized model cloud sub-ensembles. Each 
of the sub-ensembles has its own mass, heat and moisture budget. The vertical transports 
accomplished by this ensemble of model clouds is ultimately determined by the cloud base 
mass flux for each member of ensemble. In order to determine this mass flux distribu- 
tion, Arakawa and Schubert introduced the concept of quasi-equilibrium of the cloud work 
function, which leads to an integral equation making the cloud base mass flux distribution 
relating to large-scale thermodynamic processes. 

Kuo (1965, 1974) introduced cumulus parameterization for use in tropical cyclone modeling 
which were subsequently employed in large-scale numerical prediction models (e.g., Krish- 
namurti 1969; Krishnamurti et al. 1979; the NMC spectral model and the limited-area 
nonhydrostatic MM5 model). The Kuo parameterization scheme uses different mechanisms 
to describe physically the scale interaction between the cumulus ensemble and the en- 
vironment based on similar assumptions comparable to the AS parameterization scheme 
(Fraedrich 1973). The mechanism applied in the AS scheme is the vertical mass flux re- 
lating the fluxes inside and outside a convective element, i.e. , the vertical mass transport 
as the representative variable. Yet the mechanism applied in the Kuo scheme is based 
on a non-steady deep cumulus model, using the temperature difference between the cu- 
mulus cloud and the undisturbed environment and the large-scale convergence of moisture 
as indicators. The advantage of Kuo’s scheme is that, as a parameterization procedure, 
it provides immediate measures of the cumulus-scale heat and moisture fluxes in terms 
of measurable large-scale variables, without having to compute cloud dynamical processes 
(such as entrainment, detrainment and downdrafts) and cloud microphysics processes. 

Manabe et al. (1965) introduced an adjustment scheme to simulate subgrid-scale moist con- 
vection. Their adjustment includes the dry convective adjustment and the moist convective 
adjustment. The moist adjustment is performed only when the relative humidity reaches 
100% and the lapse rate exceeds the moist adiabatic lapse rate. With these artificial ad- 
justments, they introduced a simplified hydrologic cycle, which consists of the advection of 
water vapor by large-scale motion, evaporation from the surface, and precipitation process, 
to simulate the process of moist convection. 

The basic closure assumption for the AS parameterization, i.e., the cloud-work function 
quasi-equilibrium assumption, was examined by Lord and Arakawa (1980). They justified 
this assumption by considering the kinetic energy budget of a cumulus subensemble. That 
is, the generation and dissipation of kinetic energy per unit cloud-base mass flux should 
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approximately balance over time scales of the large-scale processes, and such kinetic energy 
generation (the cloud-work function) and dissipation per unit cloud-base mass flux for a 
given subensemble should not depend substantially on the large-scale conditions. They 
calculated cloud-work functions from a variety of data sets in the tropics and subtropics 
including the GATE, AMTEX, VIMHEX as well as composited typhoon data. Their results 
confirmed the correctness of the the cloud-work function quasi-equilibrium assumption. 
Lord (1982) continued to verify the AS cumulus parameterization using a semi-prognostic 
approach applied GATE Phase III data. His results show that the calculated precipitation 
agrees very well with estimates from the observed large-scale moisture budget and from 
radar observations. The calculated vertical profiles of cumulus warming and drying are 
also quite similar to the observation. His results also show that the error caused by the 
cloud-work function quasi-equilibrium assumption is generally less than 10%. In his paper, 
some additional experiments were also carried out to investigate the sensitivity of the AS 
scheme to some of the arbitrary parameters and assumptions. After these verifications and 
investigations, Lord et al. (1982) incorporated the discretized form of the AS scheme into 
the UCLA GCM. 

The AS cumulus parameterization has been widely tested and applied for various purposes. 
Ramanathan (1980) applied the AS scheme in a semi-prognostic case study of a monsoon 
depression when it was forming. He found that quasi-equilibrium held even in this disturbed 
situation, and precipitation and cumulus heating rates were reasonable. Krishnamurti et 
al. (1980) compared five cumulus parameterization schemes using the semi-prognostic ap- 
proach. The calculated rainfall rate were compared with the observed estimates from GATE 
A/B data. Moorthi and Arakawa (1985) used the AS parameterization to carry out a sys- 
tematic investigation on how cumulus heating affects baroclinic instability. Kao and Ogura 
(1987) tested the AS scheme through a semi-prognostic approach using the tropical cloud 
band data and the tropical composite easterly wave disturbance data. They found that 
cloud heating and drying effects as well as the predicted cloud population agreed rather 
well with the observations. Sud et al. (1991) tested the AS scheme using the GLA/GCM. 
They modified some parameters used in the AS scheme to improve parameterization, and 
found that the role of cumulus convection in maintaining the observed tropical rainfall and 
850 mb easterly winds was satisfactory. 

Due to the complexity of the original AS cumulus parameterization, some research efforts 
were carried out aimed at simplified parameterization schemes which would still provide 
realistic values of the thermal forcing by convection under various synoptic conditions (e.g., 
Ceselski 1974; Hack et al. 1984; Tiedtke 1989; Moorthi and Suarez 1992). Hack et al. 
(1984) used a convective flux form of the AS scheme instead of the original detrainment 
form. This flux form is more convenient to use in numerical weather prediction models. 
Tiedtke (1989) introduced a much simpler parameterization scheme and compared its results 
with the conventional convection scheme used in NWP at ECMWF. In 1992, Moorthi 
and Suarez introduced the “Relaxed Arakawa, -Schubert” (RAS) cumulus parameterization 
scheme which is a simplified AS parameterization scheme. It is very efficient, and produces 


4 



results very close to those of the standard AS implementation. The RAS scheme is used 
in the NASA GEOS-1 GCM and is described in next subsection, following closely Moorthi 
and Suarez (1992). 


2.1 Cumulus Parameterization Package of the NASA GEOS-1 GCM 


The cumulus parameterization package of the NASA GEOS-1 GCM is the RAS scheme 
proposed by Moorthi and Suarez (1992). RAS makes two major simplifications in the 
standard AS implementation: 

1) It modifies the entrainment relation to avoid the costly calculation that is required to 
find the entrainment parameter of clouds detraining at each levels. 

2) “Relaxes” the state toward equilibrium instead of requiring “quasi-equilibrium” of the 
cloud ensemble to be achieved each time the parameterization is invoked. 

The iteration method chosen in the RAS scheme considers one cloud type at a time and 
computes the cumulus mass flux that would be required to maintain the invariance of the 
work function if there were no other clouds present. It then uses the heating and drying 
effects of the given type of cloud to change the large-scale environment fields; the same 
thing is done for another cloud type based on the new modified environment fields. With 
this procedure, each step is in the direction of a single-cloud equilibrium, but in the course 
of iteration, all cloud types affect each other by modifying the environment. 

a. Cloud model 


The cloud model used by RAS is a simplified form of that in the AS scheme. The first 
major simplification introduced is to assume that the normalized mass flux for each cloud 
type is a linear function of height instead of the exponential function used in AS. Thus 

9y\{z) _ > 

dz 

where r]\(z) is the normalized mass flux for cloud type A at height z, with boundary condition 

yx(z B ) = 1 ( 2 ) 


The hydrostatic equation is used in the form 


dz 

Jp 



( 3 ) 


where P = (p/ po ) R ^ Cp , p is the pressure, R is the gas constant, c p is the specific heat at 
constant pressure, g is the gravitational acceleration, 9 is the potential temperature, and 
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Po = 1000 mb. From Eqs. (3) and (1), we have 


dy\(P) = _2 p 0X 

dP g 

and integration gives 

9x{P) = 1 + ^x1 B OdP 

9 Jp 

where Pb = P(zb), is the pressure at the cloud base height. 


( 4 ) 

( 5 ) 


As in the AS scheme, the large-scale budget of moist static energy and total water substance 
for each cloud type are 

A MPmP)] = ^jp-h(P) (6) 

and 

A {w(P) bS(P) + <S(P)]} = ^i?(P) (7) 

where h x (P), q x (P), and l x {P) are the cloud moist static energy, specific humidity, and 
liquid water mixing ratio for cloud type A at level P, respectively. h(P) and q(P) are the 
moist static energy and specific humidity in the environment, respectively. 


In Eq. (7), the precipitation term has been neglected for simplicity. It is assumed that all 
liquid water is carried to the cloud top where part is precipitated and part is evaporated, 
depending on the cloud type. Since the liquid water loading and the precipitation effects 
are excluded from Eq. (7), there is no need to specify the vertical distribution of the 
precipitation and of the liquid water. 


With these assumptions, the detrainment level of the cloud type A is the level at which the 
moist static energy within clouds equals the saturation moist static energy of the environ- 
ment. That is 

h c x (P D ) = h*(P D ) (8) 

where Pd = Pd{ A) is the detrainment level. 


Integrating Eq. (6) from Pb to Pd, we obtain 

V\(P D )h c x (P D ) - h B = -—A [ PD 9h(P) dP (9) 

9 Jp b 

Combining Eqs. (5) and (9), Eq. (8) can be solved directly for the value A corresponding 
to clouds that detrain at a given level Pp- 


MPd) 


hB ~ h*{Pp) 

C ffp*0(P)[h*(P D )-h(P)\ dP 


(10) 
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Similarly, assuming the cloud air is saturated at the level of non-buoyancy, 


q c x(P D ) = q*(P D ) 


( 11 ) 


and integrating Eq. (7) from Pjj to P#, the liquid water mixing ratio at the detrainment 
level l (P d) is calculated from 


l(P D ) = l c x (P D ) = 


Vx(Pd) 
x q(PB) + 


—A f 1 
fl Jp, 


0q(P) dP 


- q*(PD ) 


( 12 ) 


b. Cloud work function 


Following AS and neglecting the effects of water vapor and liquid water on the buoyancy, 
the cloud work function A\ is 

Ax - [ -J=n\(z)[s c x (z) - s(z)] dz (13) 

Jz B C p l (Z) 

for cloud type A, where T(z) is the temperature in the environment at height z, and s^(z) 
and s(z) are the cloud’s and the environment’s dry static energies, respectively. 


Using the hydrostatic equation (3), we have 



[s c x(P) ~ HP)] 

P 


dP 


(14) 


To obtain the cloud work function in terms of the moist static energy, let us approximate 
the static energy difference as 


s\(P) - s(P) « — * ( py [hx(P) ~ h*(P )] 


(15) 


where j(P) = ( L/c p )[dq*(P)/dT], L is the latent heat of condensation of water vapor and 
h* and q* are the saturation moist static energy energy and the saturation specific humidity 
of the environment. Finally the work function form used in the GCM is 


A\ 


f PB r,x(P) hj(P)-h*(P) 
Jp d 1 + 7 (P) P 


(16) 


c. Cumulus effects on the large-scale budgets 


The rate of change of dry and moist static energies due to cumulus convection can be written 
in the form 


/ r) q \ f) q 

[Of) = 9M c - - gLD(P)l(P)[ 1 - r(P)] (17) 
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and 


(18) 


dh\ 

~dt) 


dh 


1 = g M c ^--gD(P)(h*-h ) 


where M C (P) is the total cumulus mass flux per unit horizontal area at level P , D(P ) is 
the detrained mass per unit area and unit pressure depth, and r(P) is the fraction of the 
detrained liquid water which is precipitated. In the NASA GEOS-1 GCM (version 5.2), 
r(P) is calculated by 

r(P) = M B (\)l(P D ) (19) 


M C (P) involves contributions from all cloud types penetrating level P, i.e., 


where A (P) is given by Eq. 
mass detrainment rate 


[HP) 

M C (P) = / rj x (P)m B (X) d\ (20) 

Jo 

(10), and m B ( A) is the cloud-base mass flux per unit A. The 


D(P) = Vx (P)m B ( A) 


dX(P) 

dp 


( 21 ) 


In the NASA GEOS-1 GCM, the continuous spectrum of clouds is divided into subensembles 
of finite AA and the effects of each subensemble are considered independently. Thus for the 
*th spectral band (from A; — AA- to A;), we have 


f Ja A /-aa, Vx(P)rnB(X) dX, P > P D ( A,-), 

K(P) = Xi-AA, V\(P)m B ( A) dX, P D ( Xi) >P> P D {Xi - A A,-), 
_ 0, otherwise 


and 


D\P) 


f Vx (P)m B (X)^, P D (Xi) >P> P D (X t - AAj), 


10 , 


otherwise 


( 22 ) 


(23) 


where Pd(A 4 ) is the detrainment level of clouds with A = A i. It may be obtained by solving 
Eq. (10). 


When Xi is small, we may assume that m B (X) ~ ra#(A;), thus neglecting terms in (AA t -) 2 , 
Eq. (22) becomes 

( r] Xi (P)m B {Xi)AXi, P>P D {X t ), 

M*(P) = l Vx .(P)m B {Xi)[X(P) -X t + AA J, P D (X t ) >P> P D (X t - AA;), (24) 

( 0, otherwise 

Similarly, approximating the dX/dp in Eq. (23) by AA;[p_d(A,-) — p B {Xi — AA;)] -1 , we have 

_ f T l\, ( P) m B(X t )AX t [po(X t ) - p D (Xi - AA;)] -1 , Pd(Xi) > P > Pd (A,; - AA;), 

\ 0, otherwise 


(25) 



Then Eqs. (17) and (18) can be rewritten as 


(§0 = r *( F W A *)AA; 


and 


where 


r s(P) 


and 


r fc (P) 


= r^(P)m B (A. i )AA i 


'9V\i(P)ty 

g VXt (P)[\(P) - Xi + AX^AXi)- 1 ^ 
' +ff??A i ( j P)[p£>( A *) - Pr>(Ai - AAj)] -1 
xIa,- (Pd)L [1 - r(P D )\ , 

lo, 


P > PD{Xi), 


Pd{ Xi) >p> Pd (At - AAj), 
otherwise 


r^A.onf, 

^ Ai (P)[A(P)-A ! + AAj](AAj) -1 § 
' + 0 » 7 a,-(- p )[P£>( a *) - Pd (A* - AAj )] -1 
x[/i*(Pd) - h{Po)], 

U 


P > Pd{ X i), 


Pn(Xi) >p> pd{ Xi - AAj), 
otherwise 


(26) 

(27) 


(28) 


(29) 


Finally, the rate of change of the large-scale prognostic variables, the potential temperature 
and specific humidity due to the ith cloud subensemble can be written as 


and 


(d9\ _ m B (Xi)AX ir /D ^ 
V dt ) „ c p P sl J 


(30) 


(|) =Im B ( A s )AA s [ r ,(P)- r s ( P)] 


(31) 


d. Mass-flux kernel and cloud-base mass flux 

In AS, the change rate of cloud work function is expressed as 

dA\ _ f dA\\ ( dA x A 

dt V dt ) c V dt ) is 

where the subscripts c and Is denote the contributions from the cloud-scale and the large- 
scale processes, respectively. 



From Eq. (16) and together with some assumptions, i.e., the time variations of pressures 
at the cloud base and top equal to zero, ignoring the time dependence of 8 in Eq. (5), and 
using 


dh*{P) 

dt 


[1 + 7 (P)] 


ds(P) 

dt 


(33) 
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the change rate of cloud work function in terms of the rate of change of the large-scale 
variables, h and s, can be approximately obtained as 


dA A _ f p B dP 

~df~ Jp d P[l+l(P)] 


dh(P B ) 

dt 


[1 + 7 (P)] 


ds(P) 

dt 


+ —A 
9 



/ dh(P') 
V dt 


[1 + 7 (P)] 


ds{P) \ 
dt ) 


dP' >34) 


In AS, the cloud-scale contributions on dA\/dt is 



d\' 


(35) 


where the kernel K\ t \i represents the rate of change of the cloud work function of cloud 
type A per unit cloud-base mass flux of cloud type A'. These are not direct cloud-cloud 
interactions, but indirect effects of the various cloud types on each other through their 
environment. For quasi equilibrium to hold for the cumulus ensemble as a whole, these 
interactions must occur quickly compared to changes in the large-scale forcing. The standard 
implementation assumes that they occur instantaneously, resulting in a quasi-static balance 
between the cloud ensemble and the large-scale forcing. It is this assumption that results in 
an ill-posed problem, with the possibility that either no mass-flux distribution can produce 
an exact balance for all clouds with positive buoyancy or that (most frequently) multiple 
distributions can satisfy an “overadjustment” problem (Silva-Dias and Schubert 1977) in 
which some of the possible cloud types are overstabilized by the effects of other cloud types. 


The main assumption of RAS is that the interaction between clouds, represented by the 
off-diagonal terms of K in Eq. (35), occurs over a short but finite time and that at 
any instant the computations for each cloud and each cloud type are just based on the 
“current” environment which already has been affected by influences of some other type 
of cumulus convections. In this way, the cloud interactions are taken into account. The 
ill-posedness of the original AS implementation is thus removed by solving an initial value 
problem that selects an equilibrium distribution that depends on the time scales specified 
for the adjustment of the individual cloud types. In RAS, considering the effects of a single 
subensemble on the cloud work function, Eq. (35) reduces to 


K = 1 ( dA *A 

A ” A ’ rajg(A;)AA; \ dt ) c 

Finally, from Eqs. (34), (26) and (27), the approximate K \ it is given by 

„ f P B dP 

Ai ’ A< _ JP D P[l + l(P)] 

( r h (P B ) - [i + i(P)]r s (P) + ^a 2 [ Pb e (i v h (p ') - [i + 7 (P)]r,(P)) 

( 9 JP 



(36) 


(37) 
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The subensemble cloud-base mass flux m B (Xi)AXi is obtained by equating the large-scale 
and cloud-scale changes of A, i.e., making dA\J dt = 0, we have 





m B {Xi) > 0 , 

otherwise. 


(38) 


In the GEOS-1 GCM, the large-scale forcing (dA Xi dt)i s is calculated directly from time 
difference: 


dA Xt \ = A\ t {t + At) — A Xi (t) 
dt ) ls At 


(39) 


where A Xi (t + At) is the cloud work function calculated from the profiles of 0 and q after 
they are modified by the large-scale processes over a time interval Af. This is a good 
approximation as long as At is small. In the GCM, A Xi (t + At) is replaced by the value 
of the cloud work function after modification of the environment by large-scale effects and 
using a cloud-type dependent critical value of the work function instead of A Xi (t). Due to 
the convective parameterization is designed to nearly neutralize the instability, the critical 
value of the work function is near zero. 


In RAS, the cloud type is determined by the height of cloud-top level and the kernel is 
computed explicitly since only the diagonal elements are required. 

In addition to RAS scheme, GEOS-1 GCM employs a Kessler-type scheme for the re- 
evaporation of falling rain (Sud and Molod 1988). The scheme accounts for the rainfall 
intensity, the drop size distribution, and the temperature, pressure and relative humidity 
of the surrounding air. 

Due to the increased vertical resolution in the planetary boundary layer (PBL), the lowest 
two model layers are averaged to provide the sub-cloud layer for RAS (about 50 mb thick). 
Each time RAS is invoked (every ten simulated minutes), the possibility for shallow con- 
vection is checked for the two layers just above the cloud base. RAS also randomly chooses 
ten other cloud-top levels for the possibility of convection, from just above the cloud base 
to the model top layer. 

Supersaturation or large-scale convection is defined in the GEOS-1 GCM whenever the 
specific humidity in any grid point exceeds its supersaturation value. The large-scale pre- 
cipitation scheme rains at supersaturation, and re-evaporates during descent to partially 
saturate lower layers in a process that accounts for some simple microphysics. 

Convective and large-scale cloudiness which is used for cloud-radiative interactions are de- 
termined diagnostically as part of the cumulus and large-scale parameterizations. The 
convective and large-scale cloud fractions are combined into two separate arrays for use in 
the shortwave and longwave radiation packages. 

Supersaturation or large-scale cloudiness is defined whenever the large-scale precipitation 
scheme determines that the grid box at any level becomes supersaturated. In order to ensure 
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that at any instant the total cloud fraction is less than or equal to one, supersaturation 
clouds are only prescribed when there are no deep convective clouds. 

Since in the current moist version of GEOS-1 GCM which includes just the moist process, 
in the tangent linear model and its adjoint code the cloud formation part is not employed 
(i.e., all the related lines are commended out) although this part of code has already been 
included. 


2.2 Description of the Discretization of the Moisture Physics Parameter- 
ization 

a. Cloud model 

All clouds are assumed to have the same base. We will refer to that cloud type with its 
detrainment level in layer i, as the ith cloud type. 

The normalized mass flux for each cloud type is a linear function of height, 

Vi,k- 1/2 ~ Vi,k+i /2 = \{Zk- 1/2 - %k+ 1 / 2)5 k — i + 1, i + 2, ..., A - 1 (40) 

where r/i } k +\/2 is the cloud mass flux of the ith cloud type at level A; + 1/2 normalized by its 
value at the cloud base, X; is its entrainment rate, and Z fc+1 / 2 is the height of level A’ + 1/2. 
Equation (40) applies from the layer immediately below the detrainment layer to the layer 
immediately above the cloud base which is at K — 1/2. We assume the detrainment occurs 
at the middle of the detrainment layer and therefore that there is an additional half-layer 
at the top over which the cloud entrains, 

Vi,i ~ Vi,i+ 1/2 = -M Zi — Zi + 1 / 2 ) (41) 

The vertical coordinate is specified by the pressure at the half-integer levels (pk+ 1/21 k = 
1, 2, ..., A'). The discrete form of hydrostatic equation over the full layers is, 

Zk- 1/2 ~ Zk+ 1/2 ~ ~^k{Pk+i /2 ~ Pfc— 1 / 2 )) A; = 1,2,..., A (42) 

For the hydrostatic equation over the lower half of each layer, 

Zk - Z k+ 1 = -^0 k+1 / 2 {Pk + i - Pt). P = 1,2, ..., K (43) 

0 is the potential temperature, and P is the form as 

Pfc+ 1/2 = (Pfc+ 1 / 2 / Po) K 

and 

p _ 1 ( Pk+l/2Pk+l/2 ~ Pk-l/2Pk-\/2 

1 + K \ Pk+1/2 — Pk-1/2 



(44) 
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which is the form suggested by Phillips (1974) and used in Arakawa and Suarez (1983). 
Here k = R/c p . 

From Eqs. (40) to (42), we have 

Vi, k— 1/2 — Vi,k+ 1/2 = Pk^k\, k = i + 1, i + 2, A — 1 (46) 

where 

Pk = ~ (-Pfc+1/2 — Pk-1/2 ) (47) 

For the last half-layer up to the detrainment level, 

Vi,i = %i+i /2 + Pi°i x i (48) 

where 

Pi = ~(Pi+ 1/2 - Pi) (49) 

Summing Eq. (46) and combining with Eq. (48), we obtain the normalized mass flux at 
the detrainment level, 

i 

Vi,i = 1 + A,- ^2 PkOk (50) 

k=I\—l 

where Vi,K- 1/2 = 1 is used. 

The discrete form of the moist static energy budget of each cloud type is 

Vi,k-l/2h'~i,k-l/2 ~ Vi,k+l/2h-i } k+l/2 ~ i'Vi.k- 1/2 _ Vi,k+l/2)hk, 

k = i + 1, i + 2, ..., K — 1 (51) 

where hk is the environment moist static energy of layer k, /A fc+1/ / 2 is the cloud moist static 
energy of the *th cloud type at level k + 1/2. For the half-layer at the cloud top we have 

Vi,iK,i ~ %,f+l/2^i,i+l/2 = (%,i - Vi,i+i/ 2 )hi (52) 

Summing Eq. (51) and combining Eq. (52), we obtain the cloud-top moist static energy 
expression as 

i + 1 

Vi,iK,i = h K + 2L [(%d-l /2 _ Vi,j+ l/ 2 )^i] + (vi,i ~ Vi,i+l/ 2 )f l i (53) 

j=K — 1 

Ignoring precipitation, the cloud total water is 

Vi,k-l/2{Qi,k-l/2 + h,k- 1 / 2 ) _ Vi,k+l/2ivik+l/2 + h,k+l/2) = (Vi,k-l/2 ~ Vi,k+ 1 / 2 )%) 

k = * + 1,2 + 2, ..., A' - 1 (54) 

where g? fc+1 / 2 and l h k+ 1 /2 are the specific humidity and the liquid water mixing ratio of 
the 7th cloud type at level k + 1/2. The rhs of Eq. (54) assumes no liquid water in the 
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environment. Assuming qf i = q * , where q* is the saturation specific humidity of layer i, we 
obtain, 

l- ■ — — 

l l,l 

Vi,i 

where l h i is the liquid water mixing ratio of the detraining air for the ti ll cloud type. Since 
9 and q are known and the saturation specific humidity can be calculated and can be 
calculated if A; is known. 

Assuming at the detrainment layer 

Ki = K (56) 

where h* is the saturation moist static energy of layer i. In the AS, Eq. (56) is a polynomial 
in A ,• whose degree depends on the height of the detrainment level. From Eqs. (40), (41), 
(53) and (56) , the expression for A; is obtained as 

A h K - h* 

1 r^K-iWiW-hj) 


(57) 


i+i 

Qk + 1/2 - 1 / 2 !'/./ 1 + (vi,i - Vi,i+i/ 2 )Qi 

j=K - 1 


(55) 


b. Cloud work function 


The discrete form of the cloud work function is obtained by discretizing Eq. (16) as 


4+1 


Ai = 


J2 [ € J%j+lM h lj+l/2 - h )) 

j=K—l 


+ Mi Vi,j-l/2 ( h i,j-l /2 _ ^j)J + e iVi,i+l/2(hli +1 /2 - h*) (58) 

where A t is the cloud work function for the itli cloud type and e and /i are defined as 

c j = / ' ;+ / /2 ~ — . j = 1,2,..., K- 1 (59) 

J ^,(1 + 7 i ) ’ 

and 

M,/ = P L~ P !~y^ j = 1) 2, ..., K — 1 (60) 


PA l + 7i) ’ 

Using Eq. (51) to eliminate h c , the cloud work function can finally be written in terms of 
environmental quantities as: 


Ai — CK-ih-K — 

k 

hK + Y [(%,j— 1/2 - %,j+i/2)^i] 


*+i 

+ Yl ( e *-l “1“ Mi:) 

k=K — 1 


j=K—l 


4 + 1 


Y \ e kVi,k+l/2 + /- l kVi,k-l/2\ h * k 

k=K—l 


(61) 
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c. Cumulus effects on the large-scale budgets 


The cumulus effects on the budgets of dry and moist static energies of the environment are 
discretized as 

= ^M k -i/ 2 ( s k-i /2 — s k ) + Mk+i/ 2 { s k — s k+ 1 / 2 ) — D k lkL( 1 — r k ) j (62) 

and 

= \^k-i/2(f l k-i/2 — h k ) + Mk+\/2{hk ~ h k _ 1 - 1 / 2 ) + D k [h*k — /&&)] ( 63 ) 

where (ds k /dt) c and ( dh k /dt) c are the rate of change of dry and moist static energies of 
layer k , M k +i /2 is the cumulus mass flux at level k + 1/2, D k is the detrained mass at level 
k and A p k = Pk+ 1/2 ~ Pk-i/ 2 - The last term on the rhs of Eq. (62) represents cooling from 
the re-evaporation of liquid water detrained to the environment. Here l k is the liquid water 
mixing ratio, r k is a cloud-type-dependent precipitation fraction. 



f ds k 
V dt 


In AS implementation, M k+ \/ 2 would be the total mass flux of all cloud types penetrating 
the level k + 1/2. However in RAS, only one cloud type is considered at a time, then A7fc+i/ 2 
has the form 


M k +i/2 — 


M B {i)Vi,k+ 1 / 2 , i<k, 
0 , i > k , 


(64) 


and 


10 , 


i = k, 
i 7^ k, 


(65) 


Also l k = li t i, for i = k , and l k — 0 otherwise. Eqs. (62) and (63) can be rewritten as 


and 


where 


and 


(A) c = M fl( ,) r .W 


\ r li,k-l/2{ s k-l/2 ~ s k) + Vi,k+l/2 

r ^(A') = x(s fc - s k+1 / 2 ) - rji t ih t iL(l - r;)<5f] , k = i, i + 1, ..., K, 
0, k = 1, 2, ..., i — 1 

\vi,k-i/ 2 {hk-i /2 ~ hk) + Vi,k+ 1/2 

F/l (*0 = < x(/?. fc - h A , +1 / 2 ) + - hi) 5 f], k = i,i+ 1, ..., A, 

,0, k = 1, 2, ..., * — 1 


( 66 ) 

(67) 

( 68 ) 

(69) 
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where Sf is the Kronecker delta function and it is assumed that r/ M '_i /2 = Vi,K+ 1/2 = 0. 
Finally the rate of change of potential temperature and specific humidity due to the ith 
cloud type can be obtained as 


and 



Msji ) 

CpPk 


r ,(*), 


k = 1,2, ..., K 


(70) 


i^)=T M ^h{k)-Y s {k)l 


k = 1, 2, ..., A 


(71) 


d. Mass-flux kernel and cloud-base mass flux 


In RAS, for a single cloud type, only the diagonal element Ka is required, which is given 

by 


Ki 


1 (dAA 
M B (i ) V dt ) c 


(72) 


Using the approximate relation 


dt 


n , \ dsk 
( + ' ^ dt 


(73) 


and Eqs. (57), (61), (66), (67), 


Ki 


( £ /\ — 1 + 0)r fc (A ) — (C)?/,, ; +i /2 + d// . )( 1 + 7«)r'»(*) 


4+1 


+ i( e fc-i+lR) 

fc=A'-l l 


b /i(A ) + X! 1/2 ~ Vi,j+l/2Wh{j) 

j= K-l 


where 


+ HVi,k-i/2 - %,Jfc+i/2)rfc(*) + i?(f/,v - %i+i/2Wh(i) 


4+1 


_ X! { e kVi,k+l/2 + /lfcl?i,fc-l/2)(l + 7A:)r s (fc), 
k=K — 1 


(74) 


(%j- 1/2 - Vi,j+l/ 2 )h'j 

j=K - 1 


^ = -e t h*T] hl+ 1/2 

4 + 1 

~b yz ( £ fc-i + ^k) 

k=K — 1 

i+l 

_ [ e fc%',fc+l /2 + Mfc 7 ?*,*:— 1 / 2 ] 

k=K—\ 

i+l 

+ € ih* + i e k + /ik)hk, 

k=K — 1 


(75) 
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For simplicity, the rate of change of e and //. have been ignored. The rate of change of A; is 
included through the terms involving i ). 

With the quasi-equilibrium of cloud work function we obtain, 

M B {i) = ~ \ g K i,i t' 6 ) 

when the rhs of Eq. (76) is positive, otherwise Mb (*) = 0. The large-scale forcing of 
the Ah cloud type (Ylt) { can be calculated by Ecj. (39). Once Mg(i) is known, the 
cumulus-induced changes in 0 and q can be calculated. 

In RAS, we assume that all liquid water formed inside a cloud is carried to the top detrain- 
ment level (i.e., /j j in Eq. (55)), and a fraction of this detrained liquid water is precipitated 
and the rest is evaporated within the detrainment layer. Another assumption is the precip- 
itation simply falls to the ground without evaporation. Thus the precipitation Ri for the 
Ah cloud type can be written as 

Ri = MB{i)r t li t i (77) 

The cloud-type-dependent parameter ri is set to one for every cloud type in the NASA 
GEOS-1 GCM Version 5.2. 


3 Tangent Linear Model of the Moist Process Physics Pack- 
age 


The linearized discrete RAS parameterization equations (40) - (77) are derived as follows. 
We use {} to describe the basic state trajectory terms and ()' to denote the perturbation 
variables terms. 

For Eq. (44), the corresponding tangent linear formula is 

/ {?*+ 1/2) 4 h 

{Pk+ 1 / 2 )' = ( 78 ) 

For Eq. (45), 

( Pk )' = 1 ^ [{Pfc+1/2} - {Pfc-1/2}] 

[{Pk+ 1 / 2 } (Pk+ 1 / 2 )' ~ {pk-l/ 2 }(Pk-l/ 2 )'] 

+ 1 K [{Pfc+1/2} - {Pk- 1/2}] ({^-1/2} ~ { Pk + 1/2}) 

^{Pk-l/2}(Pk+l/2)' ~ {Pk+l/2}(Pk-l/2Y^ (79) 
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For Eq. (46) 


(%,fc-l/2) / — (Vi,k+1/2Y — {ftk}{Qk}{\)' 

k = i + 1, i + 2, /( — 1 

where 

(/?,)' =^[(n +1/2 ) , -(F fc _ 1/2 r] 


For Eqs. (48) and (49), we have 

= i 1 li,i+l/2)' + {Pi}{6i}(\y 

+ {/^* } { } (^* ) X + {6i}{K}(f3i)'i 

and 

m' = j[(Pi + i/2)'-{Pi)i 


(80) 

(81) 


(82) 

(83) 


For Eq. (50), 


(vi,iY= {a s } e [{mokY+mm'] 

k=r <- 1 

+ (Ai ) 7 Y2 {Pk}{®k} 

k=K - 1 


(84) 


The expression Eq. (55) of liquid water mixing ratio of the detraining air from the *th 
cloud type l h i can be linearized as 


(hY= -(?*)' + K0 _1 j(<7A-)' 

+ 51 {{ ( h) { r li,j-l/2) 1 - { r h,j + \/2 )'] 

i=K - i ' 

+ [{%J-l/2} - {%,j + l/2} (tfj/J 

+ (y{Vi,i} ~ {%,i+l/2}) (</i) / 

_ {ft}(%,i+l/2) , | + {%,i} 2 |{%,i+l/2}{ft} “ {^} 

51 ({%j-i/2} - H-j+i/2}) {q 3 }](vi,iY 
j=K~ 1 


(85) 
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The cloud work function expression Eq. (61) can be linearized as 
( Ai )' = {e K _i}(h K )' + {e K -\)' {hie} 

- { e i}{ 7 li,i+l/ 2 }{f l *Y 

k ( 

+ E \ [( e fc-i) ; + 

j—K—l l 


+ 

+ 

+ 

+ 

+ 


{hpc} + E ({%,j-l/ 2 } - { ? h,i+l/ 2 }) {hj} 
j=K - 1 

{efc-l} + {/^} (^/v) / 

E f ({%d-l/ 2 } - {%,i+l/2}) (fy)' 

j=K - 1 V 

(Vi,j -1/2)' ~ ( , /i,j + l/2) , j { ^ i } 


i +1 


E ■ 

k=K—l 


{ e kY{Vi,k+ 1/2} + { £ A'}(%,fc+l/2) / 


(l , fc) , { f /i,fc-i/2} + {/b’Kdi.t— 1/2/ 


TO 


{ € fc}{ ? ?j,A-+l/2} + { hi,/c— 1/2} (^fc) 


where 


[{pnaH7U)]-hw- { pj^ }) (^ 


{^•+1/2} - TO 

TOU + TO ) 2 


( 7 i) # 


and 


(M,r= - [{Pii a + { 7 ,})r (^-1/2)' 

+ r 1 -• r (, r 
tou + TO) 2 ™ 


{-Pj-l/2} 


TO 2 (l + { 7 ,» 


to 


The expressions of cumulus effects on the large-scale budgets (Eqs. (70) and 
linearized as: 


d( 9 k y 

dt 


Cp { Pk } 


(r s (k)Y 


( 86 ) 

(87) 

( 88 ) 

(71) are 
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where 

{Ap ,.} 2 {%,*— 1/2} (l s fc— 1/2} - { s t}) 

+{ r li,k+ 1/2} ({ s fc} - { s fc+l/2}) 

- {'-,}) (Ap fc )' 

+ {Ap^} iVi,k- 1/2} (( s fc-l/2) / “ (• s fc) / ) 

(r s (A))) / = < +(%,fc-i/2)' ({sfc-1/2} - {Sfc}) (91) 

+ { ? /i,fc+l/2 } (i s k - (*A'+l/2) / ) 

+ ('/i,fc+l/2) / ({ s fc} - {■‘''^+ 1 / 2 }) 

-{Vi,i}(li,iY L ( 1 - {»' })<$* 

-{Vi,iY{h,i} L ( 1 - i r i}) s i 

+L{r/ it i }{k t i } ( r t ) 'Sf , k = i, i + 1 , K, 

, 0 , k = 1 , 2 , i — 1 


20 




Finally, the cloud work function kernel and cloud-base mass flux expressions Eqs. (74) and 
(76), respectively, can be linearized as 

= (W 1 } + m)(r fc (A'))'+ ((€A--i)'+ wyr h (K)} 

- ( y { e i}{ r li,i+ I/2Y + { e iY{Vi,i+l/2} + {^}(Vi,iY + 

x (1 + {7*}){r s («)} - ( y { e i}{ r ii,i+ 1/2} + Mi 7 ?;,;}) (Tt) , {r s (*) } 

_ ({ e i}{ 7 ?i,i+i/2} + Wl 7 ?;,;}) (l + {t«}) (rs(0) / 

+ X] { ({ffc-i} + {/ 7 fc}) (r/ l (A))' 

k=K—l l 

A; r 

+ XI (i%j- 1/2} _ {^d+1/2}) (r^(j)) / 

j=K - 1 *■ 

+ (( 7 ?i,i-l/2) / - (^7 *,j-|- 1/2) {r/i(j)} 


+ ^( e fc-i) / + (/*fc)^ {T/i(A)} 

k 

+ X ({%b-l/2} “ { 7 ?i,i+l/2}) {r h(j)} 

j=K — 1 

+ {^} (i 7 ?*,*— 1/2} - { 7 /i,fc+i/2}) (IM&))' 
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+ {$} ((%,*:- 1/2)' ~ (%,fc+l/2) / ) {T h(k)} 

+ ($)' (i r )i,k- 1 / 2 } ~ {%,fc+i/ 2 }) {14(A)} j 

+ W ({»&>■} - {*7i,.-+i/2}) (14(0)' 

+ (#)' (iviA - 1+1/2}) {r fc (*)} 

+ W ((Vi,i)' ~ (%i+ 1 / 2 )') { r ^(0} 

- X { f( e fc) , { r /i,fc+l/2} + Uk}{Vi,k+l/2Y 

k=I\- > l V 

+ {/ l k}( r li,k-l/2)' + (^fc) , { ? /i,fc-l/2}^) (1 + {7 a}) {r s (At) } 

+ ({ e k}{'lli,k+ 1 / 2 } + i^k}{%k-l /2 }) (7*)'{r,(*)} 

+ ( K i e k}{Vi,k+i/A + {Mfc}{ ? ?i,fc-l/2}) (1 + {Tfc}) (r s (A)) / | (93) 


and 


where 



( 94 ) 


(*)'- 


+ 


+ 

+ 


+ 

+ 

+ 


~Ui}{ h *}(Vi,i+ 1/2)' - (Ct)'{^?}{»/M+l/ 2 } - {«.'}( fe *) , {»/M+l/2} 

i+i r 

X S ({<7—1 } + {/<fc}) 

k=K—l y 

■ a 

X [{%3- 1/2} _ {%,j+l/2} (Aj)' 

-j=A'-l 

— 1 /2) , - ( , j + 1 / 2 ) ^) {^j} 


(( c A-i)' + (fJ’k)') 


X (Kj-1/2} - {%,i+i/ 2 }) (M 1 

j=K—l J 


i + 1 

X 

k=K — 1 L 


(^{ e fc}(^7*,A'+l/2) / + ( e A:) , {%,fc+l/2} 
{/ < A-}(%,fc-l/2) / + (Ma) , {%,A-1/2}^) {^4} 
({ £ a}{^ * 7 + 1 / 2 } + {MaH^*,*- 1 / 2 }) (^4)' 


{^wy+^-mr} 
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(95) 


E 


k=K— 


1 L 


({ £ fc} + {Vk})(hiy + ((q,)' + (nkY^iK} 


3.1 Linearized Discrete Dynamical Equations 

3.2 Coding of the Tangent Linear Model 


For coding the tangent linear model, we linearize the original nonlinear forward model code 
line by line, do loop by do loop and subroutine by subroutine. This amounts to obtain the 
exact same tangent linear model as by coding directly from the original linearized model 
dynamical equations. 

The tangent linear model is the linearized nonlinear forward model in the vicinity of a basic 
state which is a model trajectory. Any original code line, we may write it as 

U = /(X) (96) 

where 

X = (aq, a; 2 , • • •, x m ) T (97) 

and U is a new derived variable related to the original control variables of the nonlinear 
forward model, i.e. , it may be one of the original control variables or an intermediate variable 
which is a function of the original control variables. Here aq, aq , • • •, x m (the components 
of the vector X) are the required variables to derive U , which may consist of either the 
original model control variables or of the intermediate variables derived from the original 
control variables; to is the number of the required variables. 


The corresponding tangent linear code assumes the form: 


SU = 


Sx i 


(df 


\dx 


+ 8 x 2 


1 X — X baste state 

df 


dx m j x =x 


basic state 


(—) 

V dxo J x— v, 

x ^ baste state 


+ ••' + 


(98) 


where X = Xbasic state means that in the expression iL, i = 1, 2, • • •, to, all the values of 
the required variables .xq, .xq, • • •, x m are chosen to have the exact same values as those of 
the basic state trajectory as in the nonlinear forward model to ensure that the basic state 
of the integration of tangent linear model is exactly the basic state of the nonlinear model 
integrating trajectory. Here. 81 _J and ^.xq, 8 x 2 , • • •, 8 x m are the corresponding perturbation 
variables of U and ,xq, .xq, •••, x m , respectively. 


In order to obtain the necessary values of X6 asJ - c state, the nonlinear model integrating 
trajectory, for the tangent linear model, we must apply the parallel method. This method 
consists of calculating in parallel the nonlinear model trajectory as the basic state Xbasic state 
and the integration of perturbation variables as well, such as <5C/, in the tangent linear model. 


23 



3.3 Notational Convention for Variables and Subroutines Used in the 
Tangent Linear Model Code 

For convenience, the same original names used in the nonlinear forward model are employed 
for the corresponding perturbation variables in the tangent linear model code. For instance, 
we use “[/” for “6U ” , “ PKHT V for U S(PKHT)”, “ AKM ” for “ 5 (ARM)”, etc.. This also 
means that the perturbation control variables in the TLM share the same common structure 
and same common block names as the full variables of the GCM itself. Thus, one needs to 
pay attention to this issue when running the TLM in conjunction with the original GCM. 

We append a “0” at the end of a variable name in the original nonlinear forward model to 
represent the corresponding basic state variable, such as “GO” for “ Ubasic state ”, “PKHT 0” 
for U (P K HT) bas i c state ” , “AKMO” for “(AKM) basic state ” , etc.. 

For naming subroutines in the tangent linear model, we simply append a “L” at the be- 
ginning of the original names of subroutines of the nonlinear forward model. To conform 
with ANSI FORTRAN 77 language, if the new name of a tangent linear subroutine exceeds 
six letters, we just retain its first six letters. For instance, the original subroutines of the 
nonlinear model “ RASG ”, “MOIST 10” and “ RNEVP ”, have corresponding names in 
the tangent linear model as “ LRASG ”, “ / MOIST” and “LRN EVP” , respectively. 


4 Adjoint Model of the Moist Process Physics Package 

4.1 Using the Adjoint Method to Calculate the Gradient of a Cost Func- 
tion 

The practical determination of the adjoint model of the moist physics package used in NASA 
GEOS-1 C-Grid GCM is the key computational method enabling us to calculate the gradient 
of a cost function with respect to initial conditions (or other control variables) for carrying 
out a 4-D variational assimilation. In 4-D variational assimilation, the cost function, which 
measures the weighted difference between observations and forecasts in an adequate norm, 
is minimized by using a large-scale unconstrained minimization method iteratively which 
requires for its implementation the gradient of the cost function with respect to the control 
variables. Finally, the optimal state defines a trajectory which passes as close as possible 
in a least-squares sense to the observations while satisfying the system of coupled partial 
differential equations of the numerical weather prediction model as strong constrains. 

Assuming that the cost function consists of a weighted least square fit of the model forecast 
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(99) 


to the observations, it has the form : 

J (x(t 0 )) = (x(tr) - X obs (t r )) T W(t r ) (x(t r ) - X o6s (t r )) 

where X(t r ) is a model state vector of size M(AK + 1) containing the values of the zonal 
wind u, the meridional wind v, the potential temperature 9 , the surface pressure P s , and 
the surface humidity q\ here M is the number of grid points at each level; K is number 
of vertical levels. t r is a given time in the assimilation window; X obs (t r ) is a vector of 
observations defined over all grid points on all levels at time t r ; W (t r ) is an N X N diagonal 
weighting matrix. From Navon et al. (1992), we have the following expression 

R T 

(VJ (X(t 0 ))) T X'(to) = E ( W (M ( X W - x °^-))) X'(t r ). (100) 

r=0 

where X'(t 0 ) is the initial perturbation, X'(t r ) is the perturbation in the forecast resulting 
from the initial perturbation, VJ (X(t 0 )) is the gradient of the cost function with respect 
to the initial conditions. 

The tangent linear model of the nonlinear forward model can be symbolically expressed as 

X'(t,,) = P r X'(t 0 ) (101) 

where P r represents the result of applying all the operator matrices in the linear model to 
obtain X'(t r ) from X'(to)- 

We define the adjoint model as 

X r (t 0 ) = PjX(t r ), r=l, (102) 

where ( ) represents an adjoint variable. After some algebra we obtain (see Navon et al. 
1992) that the expression for the gradient of the cost function with respect to the initial 
condition is 

R 

VJ (X(t 0 )) = E P^W (t r ) (x(t r ) - x° 6s (t r )) (103) 

r—O 

From this analysis, we note that the so called adjoint model operator is just the transpose 
of the tangent linear model operator. 


4.2 Coding of the Adjoint Model 

Since the adjoint model equations consist of the transpose of the linearized version of the 
nonlinear forward model, if we view the tangent linear model as the result of the multipli- 
cation of a number of operator matrices: 

P = Aj A 2 • • • Ajv, (104) 
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where each matrix A, (i = 1 , - - • , iV ) represents either a subroutine or a single DO loop, 
then the adjoint model can be viewed as being a product of adjoint subproblems 

P t = aX-!-A[. (105) 

Thus, the adjoint model is simply the complex conjugate of all the operations in the tangent 
linear model. Each DO loop or each subroutine in the tangent linear model has its adjoint 
image DO loop and subroutine, respectively. Therefore, we code the adjoint model directly 
from the discrete tangent linear model by rewriting the code of the tangent linear model 
statement by statement in the opposite direction. This simplifies not only the complexity of 
constructing the adjoint model but also avoids the inconsistency generally arising from the 
derivation of the adjoint equations in analytic form followed by the discrete approximation 
(due to non-commutativity of discretization and adjoint operations). 

4.3 Notational Convention for Variables and Subroutines Used in the 
Adjoint Model Code 

In a similar way as in the tangent linear model, we employed the same original variable 
names used in the nonlinear forward model for the corresponding adjoint variables in the 
adjoint model code. For instance, we use “A” for “A”, “PA AT” for “(PA AT)”, “ARM” 
for “(AA'M)”, etc.. As in the TLM, this convention also means that the adjoint control 
variables in the adjoint model share the same common structure and same common block 
names as the GCM itself. Thus, one needs to pay attention to it when running the adjoint 
in conjunction with the original GCM. 

We also just append a “0” at the end of a variable name (in a similar way as done pre- 
viously in the tangent linear model) to represent the corresponding basic state variable, 
such as using “AO” for “U basic stat ” , “PA'ATO” for “(PI<HT) basic state ” , “A AMO” for 
“(AA M) bas i c state ”, etc., needed in the adjoint code. 

For naming subroutines, we simply change the letter “A” at the beginning of the names of 
the tangent linear model subroutines to “A” and used them as corresponding adjoint model 
subroutine names. We also retain the adjoint subroutine names which do not exceed six 
letters to conform with ANSI FORTRAN 77 language. For instance, the original subroutines 
of the nonlinear model “ RASG ”, “MOISTIO” and “RN EVP” , have corresponding names 
in the adjoint model as “ARASG” , “AMO I ST” and “ARNEVP” , respectively. 

4.4 Verification of the Correctness of the Adjoint Model 

Integrating the nonlinear model forward in time and its adjoint backwards in time, while 
forcing the r.h.s.of the adjoint model with difference between model and observations (see 
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Eq. (103)), one can obtain values for the gradient of cost function with respect to dis- 
tributed control variables, which may consist of either the initial conditions or the initial 
conditions plus boundary conditions or model parameters. Since the moist version of NASA 
GEOS-1 C-Grid GCM consists of thousands of lines of code, any minor coding error may 
cause the final gradient of cost function with respect to the control variables to be erroneous. 
Therefore, we need to verify the correctness of the linearization and adjoint coding segment 
by segment. Each segment may consist of either a subroutine or of several DO loops. For 
a detailed derivation of the adjoint model and verification of its correctness, see Navon et 
al. (1992). 

The correctness of the adjoint of each operator was checked by applying the following 
identity (Navon et al. 1992) 


(AQ)* t (AQ) = Q* t [a* t (AQ]) , (106) 

where Q represents the input of the original code, A represents either a single DO loop or 
a subroutine. The left hand side involves only the tangent linear code, while the right hand 
side involves also adjoint code (.4*^). If equality (106) holds, the adjoint code is correct 
when compared with the TLM. In practice the identity Eq. (106) holds only up to machine 
accuracy. In our verifications of the correctness of each segment of the adjoint model and 
the whole adjoint model, the LHS and the RHS of Eq. (106) attained 13 digits of accuracy 
which is near the machine precision limit. The test were performed at NASA’s Cray C-90 
computer which has intrinsic double precision. These results show that our adjoint code 
consists of absolutely the exact adjoint operators of the TLM of the moist version of NASA 
GEOS-1 C-Grid GCM. 


A gradient check (Figs. 1, 2) was then performed to assess accuracy of the discrete adjoint 
model. This verification method is described next. First, we chose the cost function .J as in 
Eq. (99) and the N x N diagonal matrix W = diag(W u ,'W vl 'Wo,'W q ,'Wp s ), where the 
diagonal submatrices are defined as : W u = 5 X 10 — 1 1 s 2 m~ 2 1 W„ = 5 X 10 -1 I s 2 m~ 2 , 
W e = 10 _3 I A'~ 2 , W q = 5 X 10~ 3 I, Wp, s = 1 0 — 3 1 mb— 2. Then, let 

J(X + ah) = J(X) + ah T VJ(X) + 0(a 2 ), (107) 


be a Taylor expansion of the cost function. Here a is a small scalar and h is a vector of 
unit length (such as h = VJ/||VJ||). Rewriting the formula above we can define a function 
of a as 


4>(a) 


J(X + ah) - ./(X) 
ah T VJ(X) 


1 + O(a). 


(108) 
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PHI WITH RESPECT TO LOG1 0 ( ALPHA ) 



Figure 1: Variation of the <f>(a) with respect to log a (gradient check of correctness of adjoint 
model). Integration period is 6 hours and t = 6 hours model generated observations were 
used. January 1, 1985 00Z DAO’s data was used as t = 0 observations. The first guess is 
the shifted 6-hour initial condition. The time integration scheme employed is the leapfrog 
scheme. 


PHI WITH RESPECT TO L0G1 0 ( ALPHA ) 



X 


Figure 2: As in Figure 1, but for the Matsuno time integration scheme. 
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L0G1 0 ! PHI ( ALPHA )— 1 ! TO LOG1 0 ( ALPHA ) 



Figure 3: Variation of the log | <f>(a) — 1| with respect to logo-. Integration period is 6 hours 
and t = 6 hours model generated observations were used. January 1, 1985 00Z DAO’s data 
is used as t = 0 observations. The first guess is the shifted 6-hour initial condition. The 
time integration scheme employed is the leapfrog scheme. 


LOG I 0 : PHI I ALPHA )- I : TO LOG 1 0 ( ALPHA) 



Figure 4: As in Figure 3, but for the Matsuno time integration scheme. 
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For values of a which are small but not too close to the machine zero, one should expect to 
obtain values for <f>(a) which are close to unity. We obtained values for the function (f>(a) 
equal unity to a high degree of accuracy when the parameter a is varied from 10 -2 to 10 -6 . 
From the residual of (f>(a ) (Figs. 3, 4), we found that the residual tends to zero. The 
gradient check verifies that the adjoint model is correct and can be used, for example, to 
perform 4-D VDA experiments. 

Although the gradient check results are good, comparing the gradient check results with 
those for the adiabatic version of the NASA GEOS-1 GCM (Yang and Navon 1996), we 
find that including moist processes into the GCM decrease the validity of the tangent linear 
approximation. It worsens both in terms of the accuracy of the gradient values of the cost 
functional as well as in the range of the perturbations, in which satisfactory values of the 
gradient of the cost functional are maintained. The reasons are: 

1) The high nonlinearity of the original moist process package including the RAS cumulus 
parameterization and large-scale precipitation and evaporation scheme. Since the Arakawa- 
Schubert parameterization is the most complex cumulus convective parameterization scheme 
which provides the most complete physics approach and has an inherent iterative feature, 
the nonlinearity of the AS scheme as well as the RAS scheme is much stronger than that 
of some other moist parameterization schemes, such as Kuo’s scheme. 

2) In AS and RAS, the on-off discontinuous effects are more pronounced than in other 
types of moist parameterization schemes. This is due to the fact that there are more on-off 
switch processes used in Arakawa-Schubert type scheme and the iterative feature of the 
AS parameterization. Obviously these discontinuities worsen the validity of the tangent 
linear approximation and cause the values of the calculated gradient of the cost functional 
to exhibit jumps. 

Several research efforts were presented related to the serious influence of the on-off switch 
processes on the validity of the tangent linear approximation. For instance, Vukicevic 
and Errico (1993) tested the accuracy of the tangent linear model of a mesoscale model, 
compared the “true” perturbation obtained by direct nonlinear integration and concluded 
that significant errors may be expected in the regions where the moist diabatic processes 
are important for finite perturbations in the initial conditions. Recently, both Bao and 
Kuo (1995) and Xu (1996) carried out a detailed study using idealized continuous examples 
with delta function mimicking the on-off switches in physical schemes. They indicated that 
ignoring the variation of the switch point due to the perturbation in the initial conditions, 
i.e., keeping the switching point in the tangent linear model the same as in the basic state, 
could cause significant errors in the tangent linear model solution and gradient calculation. 
Thus how to deal with the on-off switches used in the moist process parameterization is a 
crucial issue. 

The essence of the influence of on/off discontinuous processes is that they may cause sudden 
jumps in the model integration trajectories, these jumps changing the model trajectories 
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and forcing the value of the cost function and its gradient to undergo a sudden change. 
This sudden change is equivalent to introducing a high nonlinearity into the variational 
assimilation system. As a consequence, it may cause either a failure or a slow-down of the 
minimization processes in 4-D Var. Simply smoothing the parameterization at the discon- 
tinuous points such as tested by (Zupanski 1993; Tsuyuki 1996a, b, c) cannot completely 
remove the influence of on/off discontinuous. Also the simple smoothing method may in- 
troduce a negative effect: that of changing the character of the original parameterization 
system. 

3) Truncation errors may worsen the validity of the tangent linear approximation. In the 
original code of the RAS scheme, there are some computational steps involving processes 
whereby a very small output results from the difference between two very large terms whose 
values are of similar magnitude, such as in calculating the rate of change of large-scale 
variables by the forcing of cumulus-scale processes. These very small time tendency terms 
result from the difference between two corresponding variable terms. Besides, a nonlinear 
term involving N dependent variables in the nonlinear forward model will create N terms in 
tangent linear computations. These newly created computation processes may increase the 
truncation error. We have checked the validity of the tangent linear approximation term 
by term and line by line in the tangent linear model, and found that the truncation error 
is a contributing factor towards a reduction in the range of validity of the tangent linear 
approximation. 
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